Contents

1 Installation

For installing cisTopic run:

devtools::install_github("aertslab/cisTopic")

In this tutorial, you will require additional packages:

install.packages('R.utils')
source("https://bioconductor.org/biocLite.R")
biocLite(c('TxDb.Hsapiens.UCSC.hg38.knownGene', 'org.Hs.eg.db'))

2 What is cisTopic?

cisTopic is an R/Bioconductor package for the simulataneous identification of cis-regulatory topics and cell states from single cell epigenomics data. cisTopic relies on an algorithm called Latent Dirichlet Allocation (LDA), a robust Bayesian method used in text mining to group documents addressing similar topics and related words into topics. Interestingly, this model has a series of assumptions that are fulfilled in single-cell epigenomics data, such as non-ordered features (‘bag of words’) and the allowance of overlapping topics (i.e. a regulatory region can be co-accessible with different other regions depending on the context, namely, the cell type or state).

cisTopic uses LDA with a collapsed Gibbs sampler (Griffiths & Steyvers, 2004), where each region in each cell is assigned to a topic based on (1) to which topic the region is assigned in other cells and (2) to which topics the regions are assigned in that cell. After a number of iterations through the data set, these assignments are used to estimate the probability of a region belonging to a cis-regulatory topic (region-topic distribution) and the contributions of a topic within each cell (topic-cell distribution). These distributions can in turn be used to cluster cells and identify cell types, and to analyse the regulatory sequences in each topic.

cisTopic consists of 4 main steps: (1) generation of a binary accessibility matrix as input for LDA; (2) LDA and model selection; (3) cell state identification using the topic-cell distributions from LDA and (4) exploration of the region-topic distributions.

Figure 1. cisTopic workflow. a. The input for cisTopic is a binary accessibility matrix. This matrix can be formed from single-cell BAM files and a set of genome-wide regulatory regions (e.g., from peak calling on the bulk or aggregate data). b. Latent Dirichlet Allocation (LDA) is applied on the binary accessibility matrix to obtain the topic-cell distributions (contributions of each topic per cell) and the region-topic distributions (contributions of each region to a topic). Note that a region can contribute to more than one topic (represented by the purple peaks). c. The topic-cell distributions are used for dimensionality reduction (e.g. PCA, tSNE, diffusion maps) and clustering to identify cell states. d. The region-topic distributions can be used to predict the regulatory code underlying the topic. For example, topics can be compared with known epigenomic signatures using a recovery curve approach; regions can be annotated and linked to genes; and, after topic binarisation, enriched motifs can be identified via RcisTarget..

If you do not want to run some of the steps in the tutorial, you can load the precomputed cisTopic object:

cisTopicObject <- readRDS('cisTopicObject_Lake.Rds')

3 Running cisTopic

3.1 Input data

Some steps in this tutorial might take some time to run, as reference we mention the running time for this dataset and settings in our system. Your actual running time will depend on your computer and dataset. In this tutorial, we will run cisTopic on 30,000 cells (with 287,000 potential regulatory regions, after removing blacklist regions) from the human brain (Lake et al., 2018).

First, load cisTopic:

suppressMessages(library(cisTopic))

The cisTopic object contains all the information and outputs from the analysis. For more information, run:

?`cisTopic-class`

In this case we will initialize the cisTopic object starting from the counts matrix. The rownames of these matrix must contain the region coordinates in position format (e.g. chr1:123456-134567). In order to save memory, we will set keepCountsMatrix = FALSE. With this option, the counts matrix will not be stored in the cisTopic object. [Reference time: 9 min]

count.matrix <- readRDS('data/counts_Lake.Rds')
cisTopicObject <- createcisTopicObject(count.matrix, project.name='cisTopic_Lake', keepCountsMatrix = FALSE)
rm(count.matrix)

By default, the slots @cell.data and @region.data will be initialized with the number of reads and accessible regions (depending on the selected threshold) per cell and region, respectively. Extra metadata can be added using the functions addCellMetadata (e.g. phenotype information) and addRegionMetadata.

cell.data <- readRDS('data/cellData_Lake.Rds')
cisTopicObject <- addCellMetadata(cisTopicObject, cell.data = cell.data)
rm(cell.data)

3.2 Building the models

The next step in the cisTopic workflow is to use Latent Dirichlet Allocation (LDA) for the modelling of cis-regulatory topics. LDA allows to derive, from the original high-dimensional and sparse data, (1) the probability distributions over the topics for each cell in the data set and (2) the probability distributions over the regions for each topic (Blei et al., 2003). These distributions indicate, respectively, how important a regulatory topic is for a cell, and how important regions are for the regulatory topic. Here, we use a collapsed Gibbs sampler (Griffiths and Steyvers, 2004), in which we assign regions to a certain topic by randomly sampling from a distribution where the probability of a region being assigned to a topic is proportional to the contributions of that region to the topic and the contributions of that topic in a cell.

To do this, runModels() builds several models (e.g. with diferent numbers of topics) using Latent Dirichlet Allocation (LDA) on the binary accessibility matrix (automatically stored in the initialized cisTopicObject). We can then select the best model using selectModel() and logLikelihoodByIter().

The main parameters for running the models (runModels) are:

  • Number of topics (topic): The number of topics are usually slightly bigger than the potential cell states in the data set. In the case of single cell epigenomics data the number of topics is low compared to other implementations (e.g. text classification). The running time will be affected by the number of topics.

  • The Dirichlet hyperparameters alpha (topic proportions) and beta (topic multinomials): Alpha affects to the topic-cell contributions; a low alpha forces to pick for each cell a few topics with significant contribution, while a high alpha allows cells to have similar, smooth topic proportions. Beta affects to the region combinations; the lower the beta, the fewer regions a topic will have; the higher the beta, the less distinct these topics will be (i.e. there will be more overlap between the topics). By default, we select alpha as 50/number of topics and beta as 0.1 (Griffiths & Steyvers, 2004).

  • Number of iterations and burnin: For recording the assignments, it is necessary that the likelihood of the model is already stabilised. cisTopic counts with the function logLikelihoodByIter to check whether this parameters should be changed. The number of iterations affect the speed of the algorithm. Note that the burnin will be substracted from the number of iterations.

NOTE: For large data sets it may not be feasible to keep all models simultaneously in memory. An alternative is to run the models and only save their likelihoods and the model with the highest likelihood (see the argument returnType in runModels). If after checking the likelihood plot another model is preferred, the function can be re-run only for that number of topics.

In this tutorial, we will test models with 2, 5, 10, 15, 20 and 25 topics. Note that running in models in parallel with such big data set can result in memory issues, and we recommend to run the models sequentially. [Reference running time for the example data set (30,000 cells; ~287,000 regions): 5 hr].

cisTopicObject <- runModels(cisTopicObject, topic=c(2, 5, 10, 15, 20, 25, 30), seed=987, nCores=1, burnin = 250, iterations = 500)

3.2.1 Selection of the best model

The log likelihood can be used to estimate the plausibility of a model parameter value, given the observed data. selectModel will select the model with the highest log likelihood (P(D|T)) at the last iteration. In order, to save memory, we can remove the binary accessibility matrix from the cisTopic object.

cisTopicObject <- selectModel(cisTopicObject, keepBinaryMatrix=FALSE)

This plot shows that a model with 25 topics is suitable. If two or more models have comparable log likelihoods, we recommend to pick the one with the lower number of topics (i.e. lower complexity). By default, this function selects the model with the highest likelihood, but the user can select a certain topic with the select parameter in this function. In cases were the topic selection is not clear, the user can rerun the models using a different seed to select the best number of topics. We will continue with the model with 25 topics.

Another way of visualizing the likelihood of the models is to plot their changes through the different iterations. It is important to check that the likelihood of the models is stabilised in the recording iterations, and the area under these curves can also be useful for model selection.

logLikelihoodByIter(cisTopicObject)

If the models are stabilized after burnin (grey line), we can conclude that the selection of the number of iterations and burnin was suitable.

3.3 Interpreting the models

3.3.1 A. Identification of cell states using the cell-cisTopic distributions

LDA returns two distributions that represent (1) the topic contributions per cell and (2) the region contribution to a topic. We can interpret these values as a dimensinality reduction method, after which the data is re-represented as a matrix with cells as columns, topics as rows and contributions as values. The recorded topic assignments to the cells (not normalised) are stored in cisTopicObject@selected.model$document_expects (see lda package).

Different methods can be used for clustering and/or visualization. cisTopic includes wrapper functions to easily run tSNE, diffussion maps, and PCA (the results are saved in the slot @dr) [Reference time: 10 min]:

cisTopicObject <- runtSNE(cisTopicObject, seed=987)
cisTopicObject <- runPCA(cisTopicObject)

Once calculations are done, cisTopic offers a unified visualization function (plotCellStates), which allows to visualize tSNE, diffussion maps, principal components and biplots (in 2/3D), colored by metadata and/or topic enrichment.

plotCellStates(cisTopicObject, method='tSNE', topic_contr='Zscore', topics='all', colorBy=c('BrainRegion', 'CombinedRegionCluster'))

We can also generate a heatmap based on the cell-cisTopic distributions.

cellTopicHeatmap(cisTopicObject, colorBy=c('BrainRegion', 'CombinedRegionCluster'))

3.3.2 B. Analysis of the regulatory topics

3.3.2.1 - Getting and visualising topics

To analyze the regions included in the cisTopics, the first step is always to derive a score that evaluates how likely is for a region to belong to a topic. getRegionsScores() calculates these scores based on the proportion of region specific assignments to a topic. These scores can be rescaled into the range [0,1], which will be useful for the binarization step (as it will force data to follow a gamma distribution shape).

cisTopicObject <- getRegionsScores(cisTopicObject, method='Zscore', scale=TRUE)

BigWig files for observing the scored regions in the genome can be generated. Note that information on the length of the chromosomes has to be provided. These files can be uploaded in IGV or UCSC for visualisation. This information can be easily found in the TxDb objects of the corresponding genomes, for example.

library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
getBigwigFiles(cisTopicObject, path='output/cisTopics_asBW', seqlengths=seqlengths(txdb))

However, many tools are limited to work with sets of regions rather than rankings of regions. Keywords, or the most contributing regions in a topic, can be used as a representative set of regions of the topic. binarizecisTopics() allows to select the top regions based on two methods:

  1. method = "Predefined": to select a predefined number of regions (determined by the cutoffs argument)

  2. method = "GammaFit" (default): to automatically select a threshold based on a fit of the scores to a gamma distribution. Note that the probability threshold must be provided by the user.

NOTE: If ‘GammaFit’ is used for the topic binarization, we recommend to scale the scores between 0 and 1 (see getRegionsScores).

cisTopicObject <- binarizecisTopics(cisTopicObject, thrP=0.99, plot=FALSE)

The regions sets selected and distributions for each cisTopic can then be analized in different ways (examples in next sections). They can also be exported to bed files to analyze with external tools:

getBedFiles(cisTopicObject, path='output/cisTopics_asBed')

3.3.2.2 - Overlap with previous epigenomic data/signatures

One way to start exploring the nature of the topics is to check their overlap (i.e. the regions included in the topics) with predefined epigenomic signatures/datasets. For example, regions from ChIP-seq can point towards enrichment of binding sites of a given TF, regions from FAIRE- or ATAC-seq to regions generally open in a given cell type or tissue, etc. In this case, we will compare our topics with the cell type specific regions identified by Lake et al. (GSE97887).

For this analysis, we need to identify which regions from the topics overlap with the signature (in this case, the whole distribution is used instead of relying only on the top regions per topic). This is achieved with an adaptation of the recovery curve approach taken by AUCell [Aibar et al., 2017]. First, epigenomic regions are intersected and mapped to regions in the dataset (by default, with at least 40% overlap). The corresponding overlapping sets (which are stored in object@signatures) are used as input, together with the region rankings per topic. Note that we can give specific labels each signature. If a region is present in the signature, the curve will increase one step in the y-axis. The area under this curve is used as a measurement for enrichment (see signaturesHeatmap).

path_to_signatures <- 'Lake_signatures/'
Lake_signatures <- paste(path_to_signatures, list.files(path_to_signatures), sep='')

labels <- sapply(strsplit(Lake_signatures, split = "Lake_signatures/GSE97887_occ1_"), "[", 2)
labels <- sapply(strsplit(labels, split = ".diffPeaks.bed"), "[", 1)

cisTopicObject <- getSignaturesRegions(cisTopicObject, Lake_signatures, labels=labels, minOverlap = 0.4)
#> The signature contains 3934 of which 3934 overlap the regions in the set.
#> The signature contains 4224 of which 4224 overlap the regions in the set.
#> The signature contains 2769 of which 2769 overlap the regions in the set.
#> The signature contains 2035 of which 2034 overlap the regions in the set.
#> The signature contains 226 of which 225 overlap the regions in the set.
#> The signature contains 12435 of which 12424 overlap the regions in the set.
#> The signature contains 4737 of which 4737 overlap the regions in the set.
#> The signature contains 28885 of which 28885 overlap the regions in the set.
#> The signature contains 10027 of which 10027 overlap the regions in the set.
#> The signature contains 7589 of which 7589 overlap the regions in the set.
#> The signature contains 18725 of which 18725 overlap the regions in the set.
#> The signature contains 8142 of which 8142 overlap the regions in the set.

We can visualize how these regions are enriched within each topic with signaturesHeatmap. With this function, we obtain a heatmap showing the row normalised AUC scores.

signaturesHeatmap(cisTopicObject)

3.3.2.3 - Annotation to genes and GO terms

Another way of gaining insight on the topics is to link the regions to genes, and to determine GO terms (or pathways or any other type of gene-set) that are enriched within them. cisTopic provides the function annotateRegions() to annotate regions to GO terms using the “TxDb” Bioconductor packages (replace ‘TxDb.Hsapiens.UCSC.hg38.knownGene’ by the appropiate organism package), and annotation databases (“OrgDb” packages).

library(org.Hs.eg.db)
cisTopicObject <- annotateRegions(cisTopicObject, txdb=TxDb.Hsapiens.UCSC.hg38.knownGene, annoDb='org.Hs.eg.db')

As we saw before, we can use the region type annotations as region sets/signatures to check whether a topic is more enriched in a certain type of region.

signaturesHeatmap(cisTopicObject, selected.signatures = "annotation") 
#> Using 4 cores.

For identifying enriched GO terms per topic, cisTopic provides a wrapper over rGREAT (Gu Z, 2017) [Reference running time: 20 minutes]. The binarised topics (i.e. sets of top regions per topic) are used in this step and results are stored in object@binarized.rGREAT. However, in this case regions are in the hg38 assembly, which is not available for rGREAT. A solution is to provide the GREAT function with a dictionary (as a GRangesList) mapping the regions in the data set with this assembly:

library(R.utils)

# url and file name for a chain file
url <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg38/liftOver/hg38ToHg19.over.chain.gz"
hg38ToHg19.chain <- "data/hg38ToHg19.over.chain"
download.file(url, destfile = paste0(hg38ToHg19.chain, ".gz"))
gunzip(paste0(hg38ToHg19.chain, ".gz"))

# Import chain file
hg38ToHg19.chain <- import.chain(hg38ToHg19.chain)

# Obtain liftOver dictionary (as list)
hg38_coord <- cisTopicObject@region.ranges
hg38_to_hg19_list <- liftOver(hg38_coord, hg38ToHg19.chain)
cisTopicObject <- GREAT(cisTopicObject, genome='hg19', liftOver=hg38_to_hg19_list, fold_enrichment=2, geneHits=1, sign=0.05, request_interval=10)

We can visualize the enrichment results within a topic:

ontologyDotPlot(cisTopicObject, top=5, topics=25, var.y='name', order.by='Binom_Adjp_BH')

We can also compare rGREAT results between topics:

ontologyDotPlot(cisTopicObject, top=5, topics=c(16,25), var.y='name', order.by='Binom_Adjp_BH')

3.3.2.4 - (Transcription factor) motif enrichment

It is also possible to identify enriched motifs within the topics and form cistromes (i.e. sets of sequences enriched for a given motif). To do this, we use RcisTarget (Aibar et al., 2017). The current version provides databases for human (hg19). You can find the region-based database at: https://resources.aertslab.org/cistarget/

For this analysis, we first need to convert the binarised cisTopic regions to the regions in the databases (“ctx regions”). This can be achieved by converting the binarised topic to a set of equivalent ctx regions (a region can map to more than one ctx region, and all regions which overlap more than the given threshold are taken). In this analysis, we also need to use the liftOver dictionary (from hg38 to hg19).

cisTopicObject <- binarizedcisTopicsToCtx(cisTopicObject, liftOver=hg38_to_hg19_list)

If needed, we can also map the regions to their corresponding ctx region (the most overlapping):

cisTopicObject <- scoredRegionsToCtx(cisTopicObject, liftOver=hg38_to_hg19_list)

We are now ready to run RcisTarget in each topic using the wrapper function topicsRcisTarget(). This function uses the binarised topic regions converted to ctx regions [Reference running time: 20 minutes]

pathToFeather <- "data/hg19-regions-1M-9species.all_regions.mc9nr.feather"
cisTopicObject <- topicsRcisTarget(cisTopicObject, genome='hg19', pathToFeather, reduced_database=FALSE, nesThreshold=3, rocthr=0.005, maxRank=20000, nCores=1)

Once RcisTarget is run, interactive motif enrichment tables can be explored (e.g. per topic):

Topic16_motif_enr <- cisTopicObject@binarized.RcisTarget[[16]]
DT::datatable(Topic16_motif_enr[,-c("enrichedRegions", "TF_lowConf"), with=FALSE], escape = FALSE, filter="top", options=list(pageLength=5))
Topic25_motif_enr <- cisTopicObject@binarized.RcisTarget[[25]]
DT::datatable(Topic25_motif_enr[,-c("enrichedRegions", "TF_lowConf"), with=FALSE], escape = FALSE, filter="top", options=list(pageLength=5))

We can also see the all the motifs in all topics:

All_motif_enr <- data.table::rbindlist(cisTopicObject@binarized.RcisTarget)
DT::datatable(All_motif_enr[,-c("enrichedRegions", "TF_lowConf"), with=FALSE], escape = FALSE, filter="top", options=list(pageLength=5))
3.3.2.4.1 - Formation of cistromes

RcisTarget results can be used to form cistromes. We define a cistrome as a set of sequences enriched for motifs linked to a certain transcription factor. In the case of cisTopic, we build topic-specific cistromes. cisTopic produces 3 different types of cistromes: ctx-regions based, original-regions based and gene based (based on region annotation). The annotation parameter decides whether only motifs linked with high confidence should be used or also motifs indirectly annotated should be considered (i.e. in this case, the *_extended cistromes will contain both annotations) [Reference running time: 8 minutes]

cisTopicObject <- getCistromes(cisTopicObject, annotation = 'Both', nCores=5)

Cistromes are useful to compare regions linked to a TF which have different spatio-temporal patterns, which may be caused i.e. by the presence of co-factors or different concentrations of the TF. Differential motif enrichment can be performed using RSAT or Homer (Medina-Rivera et al, 2015; Heinz et al, 2017); and shape features can be modelled per sequence using GBshape bigwig files (Chiu et al., 2015). These features can be used as input to Machine Learning methods (i.e. Random Forest) to determine their relevance in generating the different patterns.

Finally, you can save your cisTopic object:

saveRDS(cisTopicObject, file='cisTopicObject_Lake.Rds')

4 References

  1. Lake, B. B., Chen, S., Sos, B. C., Fan, J., Kaeser, G. E., Yung, Y. C., … & Zhang, K. (2018). Integrative single-cell analysis of transcriptional and epigenetic states in the human adult brain. Nature biotechnology, 36(1), 70.
  2. Blei, D. M., Ng, A. Y., & Jordan, M. I. (2003). Latent dirichlet allocation. Journal of machine Learning research, 3(Jan), 993-1022.
  3. Steyvers, M., & Griffiths, T. (2007). Probabilistic topic models. Handbook of latent semantic analysis, 427(7), 424-440.
  4. Aibar, S., Bravo González-Blas, C., Moerman, T., Imrichova, H., Hulselmans, G., Rambow, F., … & Aerts, S. (2017). SCENIC: single-cell regulatory network inference and clustering. Nature methods, 14(11), 1083.
  5. Medina-Rivera, A., Defrance, M., Sand, O., Herrmann, C., Castro-Mondragon, J. A., Delerce, J., … & Staines, D. M. (2015). RSAT 2015: regulatory sequence analysis tools. Nucleic acids research, 43(W1), W50-W56.
  6. Heinz, S., Benner, C., Spann, N., Bertolino, E., Lin, Y. C., Laslo, P., … & Glass, C. K. (2010). Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Molecular cell, 38(4), 576-589.
  7. Chiu, T. P., Yang, L., Zhou, T., Main, B. J., Parker, S. C., Nuzhdin, S. V., … & Rohs, R. (2014). GBshape: a genome browser database for DNA shape annotations. Nucleic acids research, 43(D1), D103-D109.

5 SessionInfo

sessionInfo()
#> R version 3.4.3 (2017-11-30)
#> Platform: x86_64-generic-linux-gnu (64-bit)
#> 
#> Matrix products: default
#> BLAS/LAPACK: /usr/lib64/libopenblas_sandybridgep-r0.2.20.so
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] doRNG_1.6.6          cisTopic_0.1.0       scatterplot3d_0.3-41
#>  [4] rtracklayer_1.38.3   AUCell_1.3.4         feather_0.3.1       
#>  [7] data.table_1.11.4    ggplot2_3.0.0        ChIPseeker_1.14.2   
#> [10] rGREAT_1.11.1        RcisTarget_1.1.2     fitdistrplus_1.0-9  
#> [13] survival_2.41-3      MASS_7.3-49          DT_0.4              
#> [16] NMF_0.23.6           bigmemory_4.5.33     Biobase_2.38.0      
#> [19] cluster_2.0.6        rngtools_1.3.1       pkgmaker_0.27.2     
#> [22] registry_0.5         fastcluster_1.1.25   destiny_2.6.2       
#> [25] Rtsne_0.13           plyr_1.8.4           doSNOW_1.0.16       
#> [28] snow_0.4-2           iterators_1.0.9      foreach_1.4.4       
#> [31] lda_1.4.2            GenomicRanges_1.30.3 GenomeInfoDb_1.14.0 
#> [34] IRanges_2.12.0       S4Vectors_0.16.0     BiocGenerics_0.24.0 
#> [37] Rsubread_1.28.1      BiocStyle_2.6.1     
#> 
#> loaded via a namespace (and not attached):
#>   [1] R.utils_2.6.0                          
#>   [2] lme4_1.1-15                            
#>   [3] RSQLite_2.0                            
#>   [4] AnnotationDbi_1.40.0                   
#>   [5] htmlwidgets_1.0                        
#>   [6] grid_3.4.3                             
#>   [7] trimcluster_0.1-2                      
#>   [8] BiocParallel_1.12.0                    
#>   [9] munsell_0.4.3                          
#>  [10] codetools_0.2-15                       
#>  [11] withr_2.1.2                            
#>  [12] colorspace_1.3-2                       
#>  [13] GOSemSim_2.4.1                         
#>  [14] knitr_1.20                             
#>  [15] rstudioapi_0.7                         
#>  [16] robustbase_0.92-8                      
#>  [17] vcd_1.4-4                              
#>  [18] VIM_4.7.0                              
#>  [19] DOSE_3.4.0                             
#>  [20] TTR_0.23-3                             
#>  [21] labeling_0.3                           
#>  [22] GenomeInfoDbData_1.0.0                 
#>  [23] bit64_0.9-7                            
#>  [24] rprojroot_1.3-2                        
#>  [25] xfun_0.1                               
#>  [26] diptest_0.75-7                         
#>  [27] R6_2.2.2                               
#>  [28] doParallel_1.0.11                      
#>  [29] RcppEigen_0.3.3.4.0                    
#>  [30] flexmix_2.3-14                         
#>  [31] bitops_1.0-6                           
#>  [32] fgsea_1.4.1                            
#>  [33] DelayedArray_0.4.1                     
#>  [34] assertthat_0.2.0                       
#>  [35] promises_1.0.1                         
#>  [36] scales_0.5.0                           
#>  [37] nnet_7.3-12                            
#>  [38] gtable_0.2.0                           
#>  [39] rlang_0.2.1                            
#>  [40] MatrixModels_0.4-1                     
#>  [41] GlobalOptions_0.0.13                   
#>  [42] splines_3.4.3                          
#>  [43] lazyeval_0.2.1                         
#>  [44] acepack_1.4.1                          
#>  [45] checkmate_1.8.5                        
#>  [46] yaml_2.1.18                            
#>  [47] reshape2_1.4.3                         
#>  [48] crosstalk_1.0.0                        
#>  [49] GenomicFeatures_1.30.3                 
#>  [50] backports_1.1.2                        
#>  [51] httpuv_1.4.3                           
#>  [52] qvalue_2.10.0                          
#>  [53] Hmisc_4.1-1                            
#>  [54] RMySQL_0.10.14                         
#>  [55] tools_3.4.3                            
#>  [56] bookdown_0.7                           
#>  [57] gridBase_0.4-7                         
#>  [58] gplots_3.0.1                           
#>  [59] RColorBrewer_1.1-2                     
#>  [60] proxy_0.4-21                           
#>  [61] Rcpp_0.12.17                           
#>  [62] doMC_1.3.5                             
#>  [63] progress_1.1.2                         
#>  [64] base64enc_0.1-3                        
#>  [65] zlibbioc_1.24.0                        
#>  [66] RCurl_1.95-4.10                        
#>  [67] prettyunits_1.0.2                      
#>  [68] rpart_4.1-13                           
#>  [69] GetoptLong_0.1.6                       
#>  [70] viridis_0.5.1                          
#>  [71] zoo_1.8-1                              
#>  [72] SummarizedExperiment_1.8.1             
#>  [73] magrittr_1.5                           
#>  [74] DO.db_2.9                              
#>  [75] SparseM_1.77                           
#>  [76] lmtest_0.9-35                          
#>  [77] mvtnorm_1.0-7                          
#>  [78] whisker_0.3-2                          
#>  [79] matrixStats_0.53.1                     
#>  [80] hms_0.4.2                              
#>  [81] mime_0.5                               
#>  [82] evaluate_0.10.1                        
#>  [83] xtable_1.8-2                           
#>  [84] smoother_1.1                           
#>  [85] pbkrtest_0.4-7                         
#>  [86] XML_3.98-1.10                          
#>  [87] mclust_5.4                             
#>  [88] gridExtra_2.3                          
#>  [89] compiler_3.4.3                         
#>  [90] biomaRt_2.34.2                         
#>  [91] tibble_1.4.2                           
#>  [92] KernSmooth_2.23-15                     
#>  [93] minqa_1.2.4                            
#>  [94] R.oo_1.21.0                            
#>  [95] htmltools_0.3.6                        
#>  [96] mgcv_1.8-23                            
#>  [97] later_0.7.2                            
#>  [98] Formula_1.2-2                          
#>  [99] DBI_0.8                                
#> [100] fpc_2.1-11                             
#> [101] boot_1.3-20                            
#> [102] Matrix_1.2-12                          
#> [103] car_2.1-6                              
#> [104] gdata_2.18.0                           
#> [105] R.methodsS3_1.7.1                      
#> [106] bindr_0.1.1                            
#> [107] igraph_1.2.1                           
#> [108] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
#> [109] pkgconfig_2.0.1                        
#> [110] bigmemory.sri_0.1.3                    
#> [111] GenomicAlignments_1.14.1               
#> [112] rvcheck_0.0.9                          
#> [113] foreign_0.8-69                         
#> [114] laeken_0.4.6                           
#> [115] sp_1.2-7                               
#> [116] annotate_1.56.1                        
#> [117] XVector_0.18.0                         
#> [118] bibtex_0.4.2                           
#> [119] stringr_1.3.1                          
#> [120] digest_0.6.15                          
#> [121] graph_1.56.0                           
#> [122] Biostrings_2.46.0                      
#> [123] rmarkdown_1.9                          
#> [124] fastmatch_1.1-0                        
#> [125] htmlTable_1.11.2                       
#> [126] dendextend_1.8.0                       
#> [127] GSEABase_1.40.1                        
#> [128] curl_3.1                               
#> [129] kernlab_0.9-25                         
#> [130] gtools_3.5.0                           
#> [131] Rsamtools_1.30.0                       
#> [132] shiny_1.1.0                            
#> [133] quantreg_5.35                          
#> [134] modeltools_0.2-21                      
#> [135] rjson_0.2.15                           
#> [136] nloptr_1.0.4                           
#> [137] jsonlite_1.5                           
#> [138] nlme_3.1-131.1                         
#> [139] bindrcpp_0.2                           
#> [140] viridisLite_0.3.0                      
#> [141] pillar_1.2.1                           
#> [142] lattice_0.20-35                        
#> [143] plotrix_3.7                            
#> [144] httr_1.3.1                             
#> [145] DEoptimR_1.0-8                         
#> [146] GO.db_3.5.0                            
#> [147] glue_1.2.0                             
#> [148] xts_0.10-1                             
#> [149] UpSetR_1.3.3                           
#> [150] prabclus_2.2-6                         
#> [151] bit_1.1-12                             
#> [152] class_7.3-14                           
#> [153] stringi_1.2.3                          
#> [154] blob_1.1.0                             
#> [155] caTools_1.17.1                         
#> [156] latticeExtra_0.6-28                    
#> [157] memoise_1.1.0                          
#> [158] dplyr_0.7.4                            
#> [159] e1071_1.6-8